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ABSTRACT 

A hot, dusty torus located around the outer edge of the broad-line region of AGNs is a fundamental ingredient in unified AGN models. 
While the existence of circumnuciear dust around AGNs at pc-scale radii is now widely accepted, questions about the origin, evolution 
and long-term stability of these dust tori remain unsettled. 

We used reverberation mapping of the hot circumnuciear dust in the Seyfert 1 galaxy NGC 4151, to monitor its temperature and rever¬ 
beration lag as a function of the varying accretion disk brightness. We carried out multiband, multiepoch photometric observations of 
the nucleus of NGC 4151 in the z, T, J, H, and K bands for 29 epochs from 2010 January to 2014 June, supported by new near-infrared 
and optical spectroscopic observations, and archived WISE data. 

We see no signatures of dust destruction due to sublimation in our data, since they show no increase in the hot dust reverberation delay 
directly correlated with substantial accretion disk flux increases in the observed period. Instead, we find that the hot dust in NGC 4151 
appears to merely heat up, and the hot dust temperature closely tracks the accretion disk luminosity variations. We find indications of 
a decreased reverberation delay within the observed period from t = 42.5 ± 4.0 days in 2010 to t = 29.6 ± 1.7 days in 2013-2014. 
Such a varying reverberation radius on longer timescales would explain the intrinsic scatter observed in the radius-luminosity relation 
of dust around AGNs. 

Our observations rule out that a second, larger dust component within a 100-light-day radius from the source contributes significantly 
to the observed near-infrared flux in this galaxy. 
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1. Introduction 


To explain the apparent dichotomy between Type 1 and Type 2 
AGNs, the unified AGN model postubtes a so-called dust torus 
iAntonuccil Il993t lUrrv & Padovanil 1 19951) around the broad¬ 
line region (BLR) of AGNs. Depending on viewing angle, 
this optically thick structure will, or will not, obscure the 
central region, resulting in different spectral energy distribu¬ 
tions (SEDs) for Type 1 a nd Type 2 AGNs. Direct and indi- 
rect observational evidence j Barvainisll 987l:lHeisler et al.lll997t 


Osterbrock & Eerlandll2006l: iJaffe et ^ 2004 ^tzuil 20061) has 

been accumulated over past decades, confirming the existence 
of the assumed dusty torus around the BLR. Nevertheless, our 
knowledge of the long-term stability of the dust around AGNs 

and of its origin remains vague. __ 

Self-consistent A.GN to rus models dKrolik & Begelm^ll988t 
ISchartmann et al.ll2010l) assume that dust is ejected into the in¬ 
terstellar medium by asymptotic giant branch (AGB) stars and 
brought from outside to the central nuclear region. Unless the 
luminosity of the central source has decreased on a timescale 
shorter than the dynamical timescale, the innermost dust for 
these models is always expected at the sublimation radius. A dif¬ 
ferent scenario, which suggests a location or formation of dust 
farther out than t he sublimation radius, has been proposed by 
lElvis et ^ (l20Q2l) . Here, the dust formation results from BLR 


clouds that are assumed to be part of an outflowing accretion 
disk (AD) wincQ. In the course of a subsequent expansion and 
cooling process, the conditions of dust condensation are fulfilled 
and dust is formed - by the AGN itself - at a distance farther out 
than the sublimation radius. 

To be able to distinguish between competing dust formation sce¬ 
narios, it seems essential to determine the stability and evolution 
of the innermost dust radius. Monitoring the hot dust tempera¬ 
ture is a further indicator of the immediate dust state, i.e. whether 
the hot dust is at the sublimation temperature. Observations, even 
if performed on one and the same AGN, NGC 4151, have so far 
given seemingly inconclusive indications: 

- V — K reverberation measurements indicate strongly varying 
time lags between 30 and 70 days over the years 2001 - 2006 
(iKoshida et al.ll20()^ . suggesting a strong variation in the lo¬ 
cation of the innermost dust over the observed period, pos¬ 
sibly by means of dust destruction and fast re-formation. A 

' It is not yet understood well dPeterson & Horn j 120061) whether the 
BLR clouds are infalling, outflowing or in rotation. Recent observa¬ 
tions indicate that the signatures of low-ionization broad emission lines 
are consistent wit h gas that is infalling or i n a rotating Keplerian disk 
dGrier et all 1201 3l : iPeterson & Hornd UOOfil) . whereas the higher ion- 
ized lines seem to trace gas that is part of an outflowing AD wind 

dCrenshaw et al.l[T9^ : lPetersonll2006l : Kollatschnv & Zetzll2013h . 
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reanalysis of these data by iHonig & Kishimot^ (1201 ih with 
a simplified clumpy torus model resulted in an apparently 
constant time lag of t = 43.8 + 8.5 days for that sa me period. 

- Interf erometric /T-band observations in 2010 dPott et al.l 
I 2 OIOI) have shown that the hot circumnuclear dust torus 
around NGC 4151 apparently does not expand owing to dust 
sublimation when the AD brightens, but found indications of 
dust survival in bright times. 

- From our 2010 z to photometry of NGC 4151, we found 
no signat ures of dust destruct ion following increased AD 
emission (ISchniille et al.l 1201^ hereafter Paper I). Instead, 
we measured a significant increase in the hot dust tempera¬ 
ture, indicating that the hot dust in NGC 4151 was apparently 
outside the current sublimation radius in that period. 

- Recently, with jif- band interferometric o bservations (Oct 
2010 - May 2012'). [iGshimoto et al.l (1201 3h have found indi¬ 
cations of an expanding /T-band dust location for NGC4151 
caused by increased AD emission and a correlation of the hot 
dust radius with the averaged AD flux over the past six years. 

While generally the dust location around AGNs is known to 
correlate with the average AD luminosity as R^ust '^Lad 
dMinezaki et alJl^Ohl : iKishimoto et al.ll201 ih . the above results 
imply that it is not a tight function of the instantaneous luminos¬ 
ity state. 

We present here the updated results from our monitoring pro¬ 
gram of NGC 4151 during 2010-2014, in which we used dust 
reverberation to study the evolution of the temperature and rever¬ 
beration lag of the hot dust around this archetypical Syl AGN. 
Compared to Paper I, we substantially extended the data set by 
using a significantly longer time series and a broader wavelength 
coverage. We also extended our reduction and analysis tech¬ 
niques by performing precise photometry with an image sub¬ 
traction method and employing a multiepoch, multiwavelength 
Markov chain Monte Carlo (MCMC) fit. NGC 4151 is part of 
our AGN hot dust reverberation campaign, in which we monitor 
roughly 25 nearby, bright, and variable Type 1 AGNs from the 
optical to the near-infrared (NIR). 

Observations, data reduction and photometry are described in 
Sect. 13 followed by an introduction of our data analysis meth¬ 
ods in Sect. [3 We present our updated results on NGC 4151 in 
Sect. |4] and discuss these in Sect. |5] A summary of our results, 
as well as the outlook for this project, is given in Sect. |3 


2. Observations and data processing 

2.1. Near-infrared data 

From 2010 January - June, we obtained NIR multiband photo¬ 
metric observation s of the Syl NGC 4 151 (see Paper I), using 
the Omega 2000 (iKovacs et al.l l2003l NIR wide-held camera 
mounted on the 3.5m telescope in Calar Alto, SpairQ. Omega 
2000 has a held of view 1514 x 1514, and its 2048 x 2048 pixel 
detector has a scale of 0"45 pixel"'. 

We continued our monitoring of NGC 4151 from 2012 February 
- 2014 June, with roughly two to four weeks of sampling (apart 
from data gaps in 2010 July - 2012 February and 2012 May - 
November). In each of the 29 epochs, we obtained broadband 
photometry of NGC 4151 and of three calibration stars in the 
same held of view in the z, Y, J, H, and K spectral bands. These 
hve hlters were chosen to differentiate between variations in the 

^ Based on observations collected at the Centro Astronomico Hispano 
Aleman (CAHA), operated jointly by the Max-Planck Institut fur As- 
tronomie and the Institute de Astrohsica de Andalucia (CSIC). 
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Fig. 1. Nuclear fluxes of NGC 4151 derived with ISIS (green) and 
GALFIT (magenta), shown here for the K band. Within the photometric 
errors, we observe excellent agreement between the two methods. The 
agreement is on a similar level for all other bands, suggesting that both 
methods manage to separate the nuclear flux of interest from the bulge 
and disk. 
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AD and hot dust emission. The weather conditions ranged from 
clear to photometric, and the seeing was 0'.'9 - 2'.'5. 

Following the data reduction with IRAF (see Paper I), we ex¬ 
tracted the nuclear fluxes of NGC 4151 at various epochs and 
filters, using two different approaches. First , we performed a 
PSF-bulge-disk decomposition with GALFIT dPeng et al.ll2002h . 
which allows for the simultaneous fitting of different galaxy 
components to the objects in an image. For NGC 4151, the fit 
model consisted of the point spread function (PSF) - constructed 
with the help of the IRAF DAOPHOT package for each frame - 
plus a Sersic profile and an exponential disk profile. This strat¬ 
egy led to small statistical errors in the fit (typically 0.01 mag 
for the PSF and Sersic components, and 0.02 mag for the disk 
component) as well as sufficiently smooth residuals for our pur¬ 
pose. The Sersic half-light radii Re, Sersic indices n, and disk 
scale lengths h for the different filters and epochs were found to 
be temporally stable. In a final fit, Rg, n, and h were fixed to their 
average values determined from the previous fits, i.e. Rg - 9'.'0, 
n - 3.0, h - 4138. These parameter value s for NGC 4151 ar e 
consistent with the literature, for example iBentz et ^ (l2009ll . 
The residuals deviate significantly from zero only in the PSF 
regime, while they are very smooth on larger scales. Substan¬ 
tial residuals in the PSF domain are likely to be caused by the 
complex PSF (due to defocusing to avoid saturation, plus see¬ 
ing), which is variable over the detector. Aperture photometry 
on these residuals yielded low fluxes, on the order of 2-3 % 
of the PSF flux of the AGN. The fluxes of the reference stars 
were determined by performing a simultaneous PSF model fit in 
GALFIT. 

Absolute flux was calibrated in the JHK bands with the 
know n fluxes of our calibr ators, as derived from the 2MASS 
PSC (ISkrutskie et al.l l2006h . For the zY photometry, we used 
a database containing model spectra of main sequence stars 
(kindly provided by R. van Boekel, priv. com.). For each cali¬ 
brator, we determined the best fit spectrum to the BVJHK fluxes 
as given by SIMBAD and calculated the zY photometry from 
that spectrum. 

As an alt ernative to GALFIT, we used the ISI S image subtraction 
package (lAlard & Luptonll998tlAi ardl2000h . which enables pre- 
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Table 1. NIR fluxes of the nucleus of NGC 4151, derived with ISIS. 


Epoch 

MJD 

z 

F. 

Y 

/ 10 

J 

H 

K 

1 

55227 

1.321 + 0.103 

1.107 + 0.074 

1.075 + 0.039 

0.970 + 0.032 

0.965 + 0.034 

2 

55252 

1.577 + 0.068 

1.324 + 0.052 

1.284 + 0.043 

1.239 + 0.045 

1.127 + 0.051 

3 

55283 

1.691 + 0.070 

1.469 + 0.050 

1.489 + 0.051 

1.499 + 0.054 

1.397 + 0.034 

4 

55311 

1.733 + 0.059 

1.550 + 0.044 

1.585 + 0.071 

1.637 + 0.052 

1.529 + 0.080 

5 

55344 

1.573 + 0.085 

1.431 + 0.058 

1.536 + 0.043 

1.723 + 0.049 

1.643 + 0.056 

6 

55373 

1.530 + 0.083 

1.390 + 0.042 

1.424 + 0.064 

1.546 + 0.050 

1.543 + 0.056 

7 

55969 

0.739 + 0.075 

0.516 + 0.056 

0.625 + 0.060 

0.712 + 0.083 

0.802 + 0.037 

8 

55998 

0.630 + 0.050 

0.483 + 0.037 

0.539 + 0.056 

0.605 + 0.040 

0.642 + 0.041 

9 

56024 

0.950 + 0.032 

0.749 + 0.039 

0.714 + 0.064 

0.653 + 0.073 

0.670 + 0.042 

10 

56255 

1.270 + 0.061 

1.038 + 0.042 

1.080 + 0.060 

1.068 + 0.040 

1.004+0.033 

11 

56264 

1.030 + 0.034 

0.881 +0.033 

0.971 + 0.031 

0.973 + 0.030 

0.969 + 0.050 

12 

56285 

1.193 + 0.082 

1.195 + 0.172 

1.089 + 0.060 

1.116 + 0.035 

1.086 + 0.039 

13 

56319 

1.169 + 0.040 

1.047 + 0.035 

1.103 + 0.033 

1.161+0.039 

1.199 + 0.037 

14 

56374 

0.953 + 0.074 

0.770 + 0.045 

0.826 + 0.064 

0.889 + 0.064 

0.915 + 0.033 

15 

56436 

1.007 + 0.030 

0.830 + 0.030 

0.901 + 0.030 

0.977 + 0.030 

0.966 + 0.030 

16 

56465 

0.721 + 0.037 

0.570 + 0.048 

0.746 + 0.039 

0.827 + 0.032 

0.913 + 0.039 

17 

56500 

0.558 + 0.035 

0.389 + 0.055 

0.526 + 0.043 

0.626 + 0.059 

0.701 + 0.039 

18 

56523 

0.613 + 0.067 

0.462 + 0.063 

0.578 + 0.057 

0.570 + 0.048 

0.648 + 0.038 

19 

56586 

0.914 + 0.039 

0.732 + 0.035 

0.791 + 0.079 

0.781 + 0.080 

0.857 + 0.048 

20 

56611 

0.917 + 0.031 

0.736 + 0.058 

0.862 + 0.032 

0.916 + 0.033 

0.909 + 0.071 

21 

56617 

0.810 + 0.034 

0.675 + 0.074 

0.792 + 0.031 

0.885 + 0.032 

0.921 + 0.040 

22 

56669 

0.500 + 0.036 

0.393 + 0.037 

0.453 + 0.055 

0.592 + 0.035 

0.729 + 0.045 

23 

56732 

0.498 + 0.039 

0.393 + 0.042 

0.378 + 0.034 

0.358 + 0.044 

0.413 + 0.051 

24 

56737 

0.559 + 0.034 

0.392 + 0.040 

0.400 + 0.039 

0.357 + 0.048 

0.413 + 0.044 

25 

56763 

0.569 + 0.040 

0.394 + 0.051 

0.453 + 0.030 

0.434 + 0.038 

0.469 + 0.034 

26 

56787 

0.505 + 0.056 

0.418 + 0.031 

0.453 + 0.040 

0.482 + 0.042 

0.520+0.037 

27 

56792 

0.501 + 0.032 

0.380 + 0.038 

0.453 + 0.035 

0.465 + 0.044 

0.533 + 0.042 

28 

56797 

0.560 + 0.043 

0.392 + 0.055 

0.464 + 0.043 

0.482 + 0.043 

0.526 + 0.042 

29 

56826 

0.551 + 0.049 

0.429 + 0.043 

0.478 + 0.047 

0.482 + 0.030 

0.532+0.030 


cise measurement of flux variati ons in variable obje c ts. Fo llow- 
ing the procedures described by IShannee & Stanelj (1201 Ih . the 
images in each ban d are first aligned using the program Sexterp 
dSiverd et al.ll2012h . Then a reference image has to be chosen or 
constructed - in our case, we chose the image with the best see¬ 
ing. A spatially variable convolution kernel K is determined that 
minimizes the discrepancy 

D = ® K]{xi,yi) - I{xi,yi)f , (1) 

i 

i.e. the kernel optimally transforms the reference image R 
to the seeing and flux level of a given individual image I. 
Each frame is then subtracted from the convolved reference 
image. Light curves, which are difference fluxes relative to the 
reference image, are then extracted from the subtracted images 
by performing aperture photometry on the residual nuclear flux. 
The error in the difference flux is estimated by measuring the 
residual flux of non-variable stars in the field (which should 
ideally be zero). Absolute calibration of the difference flux 
of each individual frame was performed with the same three 
calibration stars as in our GALFIT analysis. The absolute nuclear 
flux level in the reference image was also adopted from our 
GALFIT results. As can be seen in Fig. [T]for the K band, our 
GALFIT and ISIS fluxes agree well within the given errors. For 
targets at substantially farther distances than NGC 4151 (as are 
present in our broad AGN sample), a GALFIT decomposition 
might be more difficult owing to increased degeneracies between 
PSF and the bulge, making ISIS the method of choice for these 
objects. For this reason and overall consistency, we use the ISIS 


photometry (listed in Table [T]i throughout this paper. 

To remove emission line contributions in our zYJHK fluxes, a 
NIR spec trum of NGC 4151 was obtained with the SpeX spec¬ 
trograph (iRavner et al.ll200^ at the NASA Infrared Telescope 
Facility (IRTF) on 2010 February 26, as already described in 
Paper I. For each filter, we applied a correction factor to our data 
(ss8% in zY, ^15% in J, *2% in HK), derived from the ratio of 
the synthetic IRTF photometry to the continuum flux level at the 
respective effective wavelength. 


2.2. Mid-infrared data 

Furthermore, we used mid-infrared (MIR) images of NGC 4151 
in the WISE Wl, W2, and W3 bands, downloaded from the 
NASA / IPAC Infrared Science Archive (IRSA), and extracted 
the nuclear fluxes with GALFIT. The downloaded images in 
each band are a stack of 31 individual frames, taken on 2010 
May 31 and 2010 December 8. Although images are only 
available in stacked form, the total magnitude per band is given 
in the ALLWISE multiepoch catalog table for each of the 31 
individual observations. From a weighted magnitude difference 
between 2010 May and 2010 December (that we attribute to the 
variability of the nucleus), we derived a correction term to scale 
our nuclear fluxes determined with GALFIT to the epoch 2010 
May 31 (MJD 55347, then almost coinciding with Epoch 5 of 
our NIR observations). Our derived AGN and host fluxes are 
listed in Table |2] The total flux in each band agrees very well 
with the average total flux for NGC 4151 given in the ALLWISE 
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Table 2. WISE Wl, W2, W3 photometric fluxes of NGC 4151, derived with GALFIT. We find very good agreement with the total NGC 4151 fluxes 
as listed in the ALLWISE catalog (column 5). Our total fluxes are systematically slightly lower than the catalog fluxes, which is partly caused by 
the correction term that we applied to shift our fluxes from the mean epoch to 2010 May (see text for details). Erom 2010 May-December, the 
AGN got approximately 25% brighter. 


Filter 

AGN 

Fa / 10- 
host 

total 

total acc. ALLWISE 

Wl 

W2 

W3 

0.685 + 0.094 
0.538 ± 0.068 
0.192 + 0.035 

0.625 + 0.077 
0.301+0.101 
0.076 + 0.011 

1.311 + 0.171 
0.839 + 0.170 
0.268 + 0.047 

1.625 + 0.088 
1.022 + 0.044 
0.311 + 0.028 
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Fig. 2. Calibrated 2010 - 2013 photometry of NGC415L The optical data are plotted in orange for the epochs that coincide with our NIR epochs, 
and in black otherwise. Photometry of the NIR data was done with ISIS, and is highly consistent with our previously determined fluxes with 
GALFIT (also see Eig.[T](. 


catalog. 


2.3. Optical data 

Furthermore, optical spectra were obtained between 2012 
January and April as part of a reverberatio n mapping campaign 
described elsewhere (iDe Rosa et al.l 120141) . Observations were 
made on the 1.3-m McGraw-Hill Telescope at the MDM 
Observatory on Kitt Peak with the Boiler & Chivens CCDS 
spectrograph. A 350 lines mm“' grating yielded a dispersion 
of 1.33 A pixel”'. The entrance slit oriented north-south 
(P.A. = 0") with a projected width of 5'.'0. This configuration 
yielded spectra covering the range 4400-5850A with a spectral 
resolution of 7.9 A. We used an extraction window of 1570 in 
the cross-dispersion direction. 

Spectra were scaled to a common [O iii] /14959 flux of 
3.76 X 10”*^ ergs s”'cm”^. The starlight contribution to each 
spectrum was estimated from host-galaxy surface bright¬ 
ness models based on Hubble Space Telescope images as 
18.4 (+1.8) X 10”'^ ergs s ' cm”^ A”' (iBentz et al.ll2013h . 


The continuum fluxes at 5100 A (listed in Table O are then 
derived by averaging the spectrum over a 20 A range around the 
lowest point, which is located between the [O iii] A4959, /15007 
lines and the Fe ii blends. 


3. Methods 

3.1. Dust reverberation 

The hot dust around the central nuclear region absorbs radiation 
emitted by the AD and re-emits it in the infrared. This dust 
emission can be roughly approximated by a blackbody function. 
Owing to the absorption and re-emission of incident radiation 
by the dust, variations in the AD flux can also be detected in 
the dust emission with a characteristic lag time of t = Rdustic, 
where Rdust is the radial distance of the hot dust from the AD 
and c the speed of light. With dust reverberation, that is to 
say, by monitoring the hot dust signal as it responds to the 
AD signal, we determine the reverberation lag time and the 
temperature evolution of the hot dust around NGC 4151, to infer 
the stability of the hot dust and to conclude whether it is close 


Article number, page 4 ofllSI 



















Schniille et al.: Dust reverberation in NGC 4151 


Table 3. Optical continuum fluxes of the nucleus of NGC 4151 at 5100 A . 


MID 

Fa 

(10“'^Wm“^pm“') 

AFa 

(lO^^^Wm^^jum”') 

MJD 

Fa 

(10“'^Wm“^/im"') 

AFa 

(10^'^Wm“^/im“’) 

55932 

2.243 

0.208 

55991 

1.155 

0.204 

55933 

2.401 

0.209 

55992 

1.241 

0.204 

55934 

2.443 

0.209 

55997 

1.436 

0.205 

55935 

2.150 

0.208 

55998 

1.392 

0.205 

55936 

2.079 

0.208 

55999 

1.515 

0.205 

55937 

2.047 

0.207 

56000 

1.352 

0.205 

55938 

1.997 

0.207 

56001 

1.386 

0.205 

55940 

1.917 

0.207 

56002 

1.412 

0.205 

55941 

2.022 

0.207 

56003 

1.469 

0.205 

55945 

1.914 

0.207 

56004 

1.446 

0.205 

55946 

1.885 

0.207 

56008 

1.993 

0.207 

55947 

1.860 

0.207 

56009 

1.955 

0.207 

55948 

1.881 

0.207 

56010 

2.147 

0.208 

55949 

1.860 

0.207 

56011 

2.218 

0.208 

55950 

2.035 

0.207 

56012 

2.151 

0.208 

55952 

2.022 

0.207 

56014 

2.116 

0.208 

55954 

2.049 

0.207 

56016 

2.342 

0.209 

55957 

1.911 

0.207 

56017 

2.417 

0.209 

55958 

1.847 

0.207 

56018 

2.411 

0.209 

55959 

1.863 

0.207 

56019 

2.258 

0.209 

55960 

2.012 

0.207 

56021 

2.483 

0.209 

55961 

1.886 

0.207 

56022 

2.489 

0.209 

55962 

1.792 

0.206 

56023 

2.548 

0.210 

55963 

1.879 

0.207 

56024 

2.413 

0.190 

55964 

1.785 

0.206 

56026 

2.205 

0.189 

55968 

1.543 

0.205 

56028 

2.118 

0.189 

55969 

1.490 

0.205 

56031 

2.161 

0.189 

55978 

1.171 

0.204 

56034 

2.118 

0.189 

55979 

1.208 

0.204 

56035 

2.104 

0.189 

55981 

1.265 

0.204 

56036 

2.087 

0.189 

55982 

1.310 

0.205 

56038 

2.018 

0.188 

55983 

1.409 

0.205 

56039 

2.180 

0.189 

55984 

1.422 

0.205 

56040 

2.032 

0.188 

55985 

1.373 

0.205 

56042 

2.124 

0.189 

55987 

1.344 

0.205 

56043 

2.007 

0.188 

55988 

1.191 

0.204 

56045 

2.051 

0.188 

55989 

1.197 

0.204 

56046 

2.070 

0.189 

55990 

1.203 

0.204 

56047 

2.120 

0.189 


to sublimation. 

We note that the dust sublimation temperature depends strongly 
on the assumed dust species, composition and grain sizes, all 
of which might change significantly with radial distance from 
the source. Without detailed knowledge of these dust properties, 
it is therefore not possible to determine whether the innermost 
dust is close to sublimation from measuring the dust temperature 
in one single epoch. Monitoring the hot dust temperature for 
an enhanced period, however, allows us to draw more reliable 
conclusions about the immediate dust state. If we observe 
rising temperatures following increased AD emission, we can 
conclude that the hot dust has not yet reached its sublimation 
temperature. Clearly, the temperature evolution is a more robust 
indicator than the absolute dust temperature in one single epoch. 

To be able to measure the dust response to varying AD emis¬ 
sion, we decompose the total fluxes in each epoch and band as 
the sum of an AD and a blackbody component that we attribute 
to the dust. In contrast to the single-epoch SED fits (Eq. 13 used 
in Paper I, we perform a multiepoch, multiwavelength fit on our 


complete data set, of the form 

FAA, t, X) = Cl ■ F,(t) ■ + BB{C2, a, T{v ■ F,(t - t))) , (2) 

Here, x denotes our complete vector of model parameters. The 
first term on the righthand side of Eq. [prefers to the AD emis¬ 
sion, usually described as a power law in terms of A, with a 
power-law index a which we assume does not change signifi¬ 
cantly over the observed flux range. Eor completeness, however, 
we also test for a varying a and its influence on our results, see 
Sects. l4~2l and l5.31 Ci is a proportionality constant. In the second 
term on the righthand side of Eq.|3 BB(.) denotes the blackbody 
function 

2hc^ 1 

BB(C2,A,T)^C2-—^I^^jj^^^. (3) 

with T the blackbody temperature of the hot dust and C 2 the 
blackbody constant (consisting of a mixture of emissivity, solid 
angle, and surface filling factor). In Eq. ^ we assume that the 
AD signal has the shape of the interpolated z-band signal, which 
is justified because the z band is known to be dominated by AD 
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emission (iRiffel et al.l2009l) . Starting from an initial temperature 
To, the blackbody temperature T evolves by construction, for 
f > T, in response to the AD variations as 

dT{t)IT{t) = V ■ i ■ dL-St - T)IL,(t - r). (4) 

The variability factor v accounts for the possibility that the z- 
band variability is not com pletely reprocessed by the hot dust 
dHonig & Kishimotol201 ll) . which would lead to values of v < 1. 
Furthermore, values of v 1, might indicate that the variability 
of the heating signal is underestimated or overestimated by the 
z-band variability. The hot dust around AGN is heated mainly 
by the UV part of the AD radiation. We take the z-band flux 
as a proxy for this UV heating; however, the actual UV radia¬ 
tion is probably even more variable. Our full parameter vector 
thus comprises x = (Ci,a, C 2 , Tq, t, v). Our new multiepoch ap¬ 
proach Eq. 121141 in particular the temperature evolution given by 
Eq.m is justified by the measured hot dust temperature evolution 
obtained from preceding single-epoch SED fits of the form 


F,(T,x) = Ci 


■ + C2 • 


2hc^ 1 

^5 ^hc/AksT _ J ’ 


(5) 


i.e. AD plus blackbody contribution. These single-epoch fits 
were used in Paper I, and were also performed for each epoch 
of the complete 2010-2014 data set, prior to the multi epoch ap¬ 
proach described by Eqs.l2]l4l As can be seen in Fig. I?] the black¬ 
body temperature obtained from the single-epoch SED fits tracks 
the AD signal closely, with a stable delay of roughly 30 days, 
and the temperature variation is significant, thus justifying our 
new mode l. Eurthermore, behavior as in Eq. l4] is expected ac¬ 
cording to [Honig & Kishimotol (1201 Ih . Compared to Eq. l5j the 
data are temporally connected in our new approach. Using the 
temperature evolution described by Eq. l4]drastically reduces the 
number of parameters of the multiepoch fit, since we only need 
to estimate the initial blackbody temperature Tq instead of a tem¬ 
perature for each single epoch. This makes the algorithm more 
robust against degeneracies between T and C2. 

To check whether the distribution of the innermost hot dust is 
fairly compact or, in contrast, substantially radially extended (as 
we conjectured in Paper I), we tested for the possibility of a fur¬ 
ther blackbody component contributing to the hot dust signal. 
Therefore, in addition to the IBB-model of Eq. |2l we fitted a 
2BB-model of similar type to our data: 


F,i{A,t,x) = Ci-F,{t)-A-‘^ + BBliC 2 ,A,Ti{vi-F,{t-Ti))) 

+BB2(C2, a, T2(v 2 ■ F,(t - T2))) (6) 


with the temperature evolution of the second blackbody 
given by an expression as in Eq. |4l though with differ¬ 
ent starting temperature, time lag, and variability factor. 
This 2BB-model has four additional parameters, and x = 
(Cl, O', C2, To,i, Ti, vi, C3, To,2, T2, V2)- 


3.2. Interpolation of the AD signal 

To calculate the blackbody temperature T described by Eq.|4]of 
our new approach for any arbitrary value of t in the range of 
the data, we first need to interpolate our input accretion disk sig¬ 
nal, i.e. the z-band signal. This is done with the method for in¬ 
terpolation, real ization, and reconstruct ion of noisy, irregularl y 
sampled data bv iRvbicki & PressI (Il992h and iPress et al.l (i l992h . 
which was pre viously used for re verberation studies of the BLR, 
for example bv IZu et ^ (1201 ll) or lHernitschek et ^ (l2014li . The 



T / yrs 

Fig. 3. Estimated (v,j)'^^ (red dots) from the z band data versus T,y, 
shown on logarithmic scale. The three different structure function mod¬ 
els that we fit to these data (see text of Sect. 13.2t are represented by the 
green, magenta, and blue lines. 


AD signal can be described as a stochastic process character¬ 
ized by two structure function parameters, which are estimated 
from the data. With this we can predict the signal for unmeasured 
times. 

Our M - 29 measurements y = (y,), / = 1,..., M of the z-band 
flux can be written as 

y = s H- Ey H-n, (7) 

where s represents the intrinsic variability signal, E is the vec¬ 
tor (1,1,1..., 1)^, y is an appropriate estimator of the mean of 
the data, and n represents the measurement noise. The intrinsic 
signal s can be described as a stochastic process with a certain 
correlation structure, and we can calculate an estimate s* of the 
signal at any (measured or unmeasured) point using the already 
obtained M measurements: 

s. = d^(y-Ey)+ . (8) 

Here, d, are linear coefficients that depend on the particular point 
to be estimated, and x, is the discrepancy between the estimated 
and the true value. Minimizing the discrepancy with respect to 
the linear coefficients then yields the solution for the coefficients, 
and the least variance estimate of our signal 

T. =Sr[S+N]-'(y-Ey)+y. (9) 

Here, S is the covariance matrix of the data, S, is the covariance 
vector between the measured data points and the new data point, 
and N is the noise covariance matrix. In Eqs.[8]and|9] the mean^ 
is subtracted from the data before determining the coefficients 0 
and is added back to the estimate afterward. The variance of the 
true value i, about the best estimate T, is given by 

<(s, - s;)2> = <s2) - sHs + N]-'s.. (10) 

To apply Eq. |9] we need to estimate the covariance structure 
of the underlying stochastic process from our data. Assuming 

^ The appropriate mean is calculated by 

. _ E^[S-tN]-‘y 
E^[S + N]-‘E’ 

and corresponds to the value y that minimizes 
;t^" = (y-Ey)^[S + Nr‘(y-Ey) 


when subtracted from the data. 
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stationarity, so that S ij = S (f,- - tj) = S (r), the covariance matrix 
can be replaced by the population mean square and the structure 
function V(t): 

( 11 ) 


V(t) = + (12) 

Following the method described by iPress et ^ (Il992ll . we cal¬ 
culate a lag Tij - \ti -tj\ and an estimate of the structure function 
t'lj - (yi-yj)^ “ j pair (i, j), where n,- is the mea¬ 

surement error of data point y, . The M(M - l)/2 = 406 pairs are 
binned into 20 bins, equally spaced in r. In Fig. [2 is plotted 

against r/j on a logarithmic scale. We fit a straight line 

log( = log(A) -H |log(T) (13) 

to these data, representing a power-law model of the form 


where A is the average flux variability on a one-year timescale, 
and y is the gradient of this variability. For the structure function 
parameters, we obtain best fit values of A = 0.011 + 0.001 Jy 
(resp. A = 0.423 + 0.019W and y = 0.881 + 0.067. 

With these values, we calculate the covariance matrix and 
vector, and interpolate ourz-band light curve according to Eq.|2 
The interpolated z-band light curve is shown in Fig.|4] Since our 
measurement errors are generally small (at the 5% level), the 
interpolated light curve runs almost perfectly through our data 
points. In Fig. [3j a potential substructure seems to be evident 
around r = 2.5 yr. However, this apparent feature might be 
merely an artifact, caused by insufficient sampling in that range 
of T. Indeed, the bins around t = 2.5 yrs contain by far the 
least number of data points, obviously owing to the large data 
gap between June 2010 and February 2012. An apparent feature 
seen in Fig. [3] seems to be common in the empirical s t ructur e 
functions of si n gle tar g ets (see e.g. figure s in Press et aT] dl992h : 
ISchmidt et al.l OOlOh : iMorganson et al.l (12014 )'). while it is 
averaged out i n the observed struct ure function of samples of 
AGNs (see e.g. lSchmidt et alJ(l2010h '). 

Our derived structure function parameters agree with values 
found in the literature. While for quasars, typical values of 
the power-law slop e are y a; 0.3 - 0.4 (iBauer et al.l 120091 : 
ISchmidt et al.ll201^ . the structure functi on of Seyferts is gener¬ 
ally found to be steeper (lHawkinsll2002h . except for very short 
timescales, wh ich are dominated b y noise. Indeed, specifically 
for NGC 4151. I^ernv et al] (l2003h report slopes of a; 0.65 - 1.0 
for the V band structure function on different timescales and in 
different epochs. 


An alternative approach of the structure function to describe 
AGN variability is the damped random walk (DRW) model: 

y(T) = 0-^(1- exp(-T/T^)), (15) 

where cr^ is the long-term variance of the process, and the 
damping timescale. While for a substantial t range the variabil¬ 
ity increases with increasing time lag following a power law 
V{t) cx t^, the DRW model considers that there is a plateau 
where the variability saturates, at time lags that are longer than 
the longest correlation timescale of the stochastic process. When 
we fit this model to our observed structure function for the NGC 



MJD 


Fig- 4- z-band ligh t c urve, interpolated with the method of 
iR^^icki & PressI ( 1199^ and iPress et all ( 119921) . We interpolated using 
the three models described in Sect. 13.21 For clarity, only the interpo¬ 
lation resulting from the power-law structure function is shown here. 
This interpolated signal is the least-variance estimate of the stochastic 
process and therefore very smooth. The error of this prediction is indi¬ 
cated by the shaded gray region. Single realizations of the process have 
much more structure on short timescales and will also make excursions 
outside of the estimated error regions. 


4151 z band signal, we obtain values cr - 0.148 + 0.505 Jy and 
Td = 230.0 + 576.1 yrs. Clearly, the fitted damping tim escale, 
typically between 0.1 and 3 years dMacLeod et akll^lOl) . seems 
physically unreasonable. When looking at the observed structure 
function data plotted in Fig. |2 we can see that our time series is 
obviously not long enough to sufficiently constrain the damp¬ 
ing timescale for NGC 4151, as also evident from its large error. 
Obviously, no plateau is reached, but the variability increases 
througho ut the observed time range. This is consistent with the 
results of ICzernv et al.l(l2003h . who find an unusually long damp¬ 
ing timescale of roughly ten years for this object. We performed 
the fit again with fixed Td - IQ yrs, leading to cr = 0.032 + 0.002 
Jy. 

The obtained structure functions for the power law and the DRW 
model are plotted in Fig. |3 To test the influence of the particular 
choice of structure function on the robustness of our results, we 
performed all our following analyses (i.e., the interpolation of 
the AD signal as well as the various fits presented in Sects. 14.21 
andQ with all three models shown in Fig. [3] 

3.3. Model fitting using a DE-MC aigorithm 
Differential evolution Markov chain 

Compared to Paper I, in which we used the routine mpfit to fit 
the model of interest to our data, we now infer the best model pa¬ 
rameters using a differential evoluti on Markov chain (DE-MC) 
algorithm, which was proposed by iTer Braafl (l2006h . DE-MC 
is a population MCMC algorithm, with multiple chains running 
in parallel. Compared to random walk metropolis (RWM) type 
MCMC algorithms, the main advantage of DE-MC is that the 
problem of choosing a jumping distribution, hence an appro¬ 
priate scale and direction of the jumps, is solved by generating 
jumps as a fixed multiple of two randomly selected parameter 
vectors Xri and Xs 2 that are currently in the population: 

Xp = X; -H 7 ■ (x«i - Xr 2 ) + e (16) 

where y > 0, and the vector e of small random numbers is added 
to ensure that the whole parameter space can be reached. These 
random numbers are drawn from a symmetric (e.g. Gaussian) 
distribution that has small variance compared to the target 
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variance. In the proposal scheme of Eq. [16] the chains learn 
from one another, hence the name “differential evolution”. 




2 1.3 - — Chain 5 

\ — Population mean 

= 1 . 2 - - 

I— 

1.1 _^^ 



^ 0.0100 
o 

O 

^ 0.0010 
> 0.0001 



Fig. 5. Evolution of the parameter Toj (for one chain and the popula¬ 
tion mean) of the IBB model, for one exemplary run. Also shown are 
the evolution of the pooled within chain variance W (Eq. IA.4t . maxi¬ 
mum variance estimate V (see Eq. IA.lt and the resulting convergence 
parameter, i.e., the PSRE R (Ea. IA.lt . Convergence was diagnosed at 
iteration step n = 13599, where R < 1.1 was reached for all parameters. 

Implementation of the algorithm 


We ran the DE-MC algorithm for the six-parameter IBB 
model of Eq.|2| with 15 simultaneous Markov chains and 25 si¬ 
multaneous chains for the case of the ten-parameter 2BB model 
of Eq.| 6 | We found y - 2.38/ to be a good value for achiev¬ 
ing the desired acceptance rate of roughly 25%. In each iteration 
step and for each chain, the algorithm evaluates - for the current 
parameter vector x, and for the proposal vector Xp - the posterior 
probability density function (pdf) 

p(x|y) oc p{x) ■ p(y|x) 

CX p(x).exp(-^| 


which is proportional to the prior pdf times the likelihood func¬ 
tion according to Bayes’ theorem. Here, x represents the param- 


Epoch3, T=1471K Epochs, T=1497K 






X //um X / fJ-m 

Fig. 6. Our data and the IBB model fit for the epochs 3,5,9, and 10. 
Overplotted is the AD model, the IBB model, and the sum of both com¬ 
ponents. From the data, we observe a clear rise in the blackbody flux, 
and a clear shift of the blackbody emission peak to shorter wavelengths 
until epoch 5. This emerging bump is confirmed by the blackbody tem¬ 
perature maximum in epoch 5 of our single-epoch fit (see Fig.jT}. 

eter vector, and y are the data. The notations and respec¬ 
tively refer to the values of the model and the data for the kth 
of the = 211 NIR plus optical data points (the WISE data 
were not included in the model fit), and cric is the related photo¬ 
metric error. To not overweight the low flux epochs around MJD 
56000, where all of the 5100 A measurements are assembled, we 
decided to use a weighted= xl/m + Asioo/^- Thus, we obtain 
an effective number of A = 160 data points. 

Eor each parameter x, we chose a uniform prior pdf 

f 7^, a < X < b 

P(x) = ] (18) 

(^ 0 , X < a \/ X > b 

within the limits a,b as given in Table 01 in order to exclude 
unphysical parameter ranges. The posterior pdf is then directly 
proportional to the likelihood, and the maximum of the posterior 
occurs when is smallest. The upper limit for Ci is given by 
the fact that for the z band (T 0.9pm), Ci ■ 2“ must not exceed 
1. In the 2BB model, the maximum temperature for the second 
blackbody component is given by T^ax - T\ ■ yfrijr^, corre¬ 
sponding to the case that all radiation from the AD can reach the 
second blackbody, i.e. the covering factor of the first blackbody 
is negligible. We chose T 2 = 100 days as an upper limit for the 
reverberation lag of the second blackbody component, because 
at that distance (and resulting low temperatures), its contribution 
to the NIR fluxes would become negligible, unless the black¬ 
body constant C 3 would be several orders of magnitude higher 
than C 2 . We further set p(t 2 ) - 0 for T 2 < ti . 

As initial distribution for the parameters, we took a uniform dis¬ 
tribution within the limits of our priors, which is thus sufficiently 
overdispersed. 

Eor monitoring convergence of the D E-MC algorithm, we used 
the convergence criteria prop osed by iGelman & RubinI (Il992h 
and iBrooks & GelmanI dl998l) . i.e., the R value, and visual in¬ 
spection of W and V, as explained in Sect. [^ The IBB model 
converged within 10000-15000 iteratio ns on average. We con¬ 
sidered R < 1.1 as close enough to 1. iTer BraakI (120061) uses 
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Table 4. Upper and lower limits a,b for uniform pdfs of the model parameters, for the IBB model (with constant and varying a] and the 2BB 
model. 


Parameter 

Unit 

IBB 

model 

IBB 

model (var. a) 

2BB 

model 



a 

b 

a 

b 

a 

b 

Cl 

- 

0.0 

0.85 

0.0 

0.85 

0.0 

0.85 

a 

- 

1.5 

2.5 

1.5 

2.5 

1.5 

2.5 

ai 

- 

- 

- 

0.0 

1.0 

- 

- 

C 2 

lO-'^ster 

1.0 

30.0 

1.0 

30.0 

1.0 

20.0 

To.i 

lOOOK 

1.0 

2.0 

1.0 

2.0 

1.0 

2.0 

Tl 

days 

0.0 

80.0 

0.0 

80.0 

0.0 

40.0 

Vl 

- 

0.5 

2.0 

0.5 

2.0 

0.5 

2.0 

C3 

lO-'^ster 

_ 

_ 

_ 

_ 

10.0 

100.0 

To,2 

lOOOK 


_ 

_ 

_ 

0.5 

1.0 

T2 

days 

- 

- 

- 

- 

30.0 

100.0 

V2 

- 

- 

- 

- 

- 

0.5 

2.0 


Table 5. Global mean and errors for the parameters of our various models, and reduced of the fit. The errors are given by Ax = xfV for each 
parameter. All results shown here are derived from using the power-law structure function. 


Parameter 

Unit 

IBB 

IBB (var. a) 

IBB (ep. 1-9) 

IBB (ep. 7-29) 

2BB 

Cl 

- 

0.78 + 0.01 

0.78 + 0.01 

0.78 + 0.01 

0.78 + 0.01 

0.77 + 0.01 

a 


1.63 + 0.04 

1.56 + 0.07 

1.65 + 0.07 

1.63 + 0.04 

1.65 + 0.04 

ai 

- 

- 

0.10 + 0.07 

- 

- 

- 

C2 

lO-'^ster 

3.82 + 0.19 

3.88 + 0.20 

4.63 + 0.33 

3.52 + 0.22 

3.16 + 0.34 

To.i 

lOOOK 

1.436 + 0.015 

1.436 + 0.015 

1.392 + 0.020 

1.455 + 0.018 

1.479 + 0.027 

Tl 

days 

31.0+ 1.6 

30.8+ 1.5 

42.5 + 4.0 

29.6+ 1.7 

29.3+ 1.9 

Vl 

- 

0.81+0.02 

0.82 + 0.02 

0.96 + 0.04 

0.81+0.03 

0.80 + 0.02 

C3 

lO-'^ster 

- 

- 

- 

- 

56.47 + 25.65 

To,2 

lOOOK 





0.698 + 0.064 

T2 

days 

- 

- 

- 

- 

67.2 + 8.3 

V2 

- 

- 

- 

- 

- 

1.26 + 0.20 


- 

1.87 

1.87 

1.13 

1.76 

1.89 


R < 1.2. Figure |5] shows the evolution of the parameter Tqj 
(which took longest to converge) for one exemplary run, as well 
as the evolution of R, V and W for this parameter. Visual in¬ 
spection shows that W and V are not evolving anymore and con¬ 
vergence was diagnosed correctly. The alternative ten-parameter 
2BB model took 50000-100000 iterations on average to con¬ 
verge. 

4. Results 

4.1. Constant power-law index 

All results presented in this section were obtained using the 
power-law structure function model and a constant AD power- 
law index a. For results under the use of different structure 
function and a models, see Sects. I4.2I and I5.3I The temporal 
evolution of the photometry of the nucleus of NGC 4151 is 
shown in Fig. |2] The flux variations in the other bands show 
a time lag behind those of the AD-dominated z band - except 
for the 5100 A fluxes, which seem to be concurrent with the 
z-band fluxes - and this lag increases with wavelength. This 
band-dependent behavior can be explained by the reverberation 
delay of the hot dust and by the increasing dust contribution and 
decreasing AD contribution for longer wavelengths. 

Figure [T] shows a single-epoch decomposition of the fluxes 
into AD and hot dust according to Eq. |5] The hot dust tem¬ 
perature closely follows the overplotted AD z-band flux with 
a delay of 30 days, which justifies the model described by 


Eqs. IHH] The minimum that is seen in the temperature around 
MJD 56700 (and also in the corresponding JHK flux minima) 
is not covered within the sampling of the z-band data (see Eig.[ 8 ] 
and Eig.|2]i, therefore we observe a substantial discrepancy here. 
Comparing epochs 1-9 with epochs 10-29, one can see that the 
ratio Fad.zIT is higher in the first epochs, already indicating 
our later finding that in these early epochs the reverberation lag 
might be larger than in the later epochs, so that the AD flux is 
more diluted and the dust heated less in epochs 1-9. 

The results of our multiepoch, multiwavelength fit is shown 
in a temporal plot in Eig. [8] for the IBB model and in Eig. 0 
for the 2BB model. Eor the IBB model, the best-fit lag is 
T - 31.0 + 1.6 days, and the best-fit value for the power-law 
slope is a = 1.63 + 0.04 (cf. Table |5] for values of the other 
parameters), nicely matching the Fy oc law expected from a 
standard Shakura-Sunyaev accretion disk. The initial blackbody 
temperature has a best-fit value of Tq = 1436 + 15 K, and then 
evolves (by construction) as in Eq. |4l Epochs 1 and 2 are not 
included in the evaluation of the fit, as in our model, T evolves 
only for f > r. We find;ifj^^ = 1.87 for this fit, indicating that we 
have probably slightly underestimated our measurement errors 
or that the model is too simple. 

Eor the 2BB model, we And best-fit values of ti = 29.3 +1.9 
days, To,! = 1479 ± 27 K, T 2 = 67.2 + 8.3 days, To ,2 = 698 ± 64 
K,a- 1.65 + 0.04, and a reducedvalue - 1-89. Here, 

epochs 1-3 are not included into the fit, because T 2 evolves only 
for f > + 2 . The best-fit values for the complete set of parameters 
are given in Table |5] Eor both the IBB and the 2BB model, the 
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Fig. 7. Hot dust temperature evolution derived from the single-epoch decomposition. The hot-dust temperature changes follow the relative AD 
z-band flux changes with a delay of roughly 30 days. 




MJD 






MJD MJD 


Fig. 8. Results of our IBB model fit. Plotted are the K, H, J, Y, z, and 5100 A data (green) over time. The bigger dots with errors bars in the 5100 
A panel mark the data points that coincide with our NIR epochs. The red line represents the blackbody contribution for each band, blue the AD 
contribution, and the black line is the sum of both. Epochs 1 and 2 are not included in the evaluation of the fit (as t < t). 


data can be fit well with a stable time lag. 

In Fig.|6] we show the results from our IBB fit in the spec¬ 
tral domain, for the epochs 3, 5, 9, and 10. It is clear how the 
observed SED changes from epochs 3 to 5, and we can see a hot 
dust bump emerging, with the peak of the blackbody emission 
shifted to lower wavelengths (A » \.9pm), thus indicating our 
detected temperature increase. Though degeneracies between 


the blackbody constant and the blackbody temperature are non- 
negligible (in the single-epoch fits, where the temperature can 
evolve freely, the correlation coefficient is as high as -0.75), 
visual inspection of the data and the systematic changes in 
the NIR color, correlating with the delayed AD brightness, 
underline the actual temperature increase. The hot dust peak is 
shifted to much longer wavelengths (T a; 2.5pm) in epoch 9, 
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Fig. 9. Results of our 2BB model fit. Plotted are the K, H, J, Y, z, and 5100 A data (green) over time. The bigger dots with error bars in the 5100 
A panel denote the data points that coincide with our NIR epochs. The red line represents the contribution of the inner blackbody component (at 
29 light-day distance) for each band, the brown line is the second blackbody (at 67 light-days), blue the AD contribution, and the black line is the 
sum of all three components. Epochs 1, 2, and 3 are not included in the evaluation of the fit (as t < T 2 ). The resulting 2BB model fit only differs 
significantly from the IBB model fit in the high-flux epochs 4-6, in the K band, as indicated by a red circle. 


and until epoch 10, we again observe rising temperatures. 


In Fig. [TOl we show the marginalized posterior probability 
distributions for the parameters C 2 , Tq, t, and a. As can be seen, 
C 2 and To are strongly anti-correlated. However, any increase 
or decrease in C 2 , hence decrease or increase in Tq, will only 
shift the resulting temperature curve in the vertical direction, 
owing to our approach Eq. |4l The significance of the tempera¬ 
ture variations justifying this approach has already been shown 
with our single-epoch fits (Fig.|2l) and can also be seen from the 
SEDs shown in Eig.|6l The reverberation delay r does not seem 
to show significant correlations with any other parameter, while 
a is slightly anti-correlated resp. correlated with C 2 resp. Tq. 
Interestingly, we observe a multimodal pdf of the reverberation 
lag T. This is on the one hand caused by slightly different re¬ 
verberation delays in 2010 and 2013-2014 (see Sect. 15.3|) . but 
mainly by a strong bimodality in r in the 2012-2014 part of the 
data set (epochs 7-29, see Fig.ITTI). Because that second period 
dominates the fit due to a higher amount of data points and lower 
photometric errors, the bimodal pdf is also visible in the fit of the 
complete data set. This bimodality is discussed in more detail in 
Sect. 153] 


4.2. Further structure function models and the time-variable 
power-law index 

As described in Sect. 13.21 we tested for the robustness of all our 
performed multiepoch, multiwavelength fits using three different 
structure function models. It turned out that the results are highly 
stable under the exchange of the structure function model or its 
particular parameters. In Table |5] the results are only listed for 
the power-law model, while Table |6] shows the influence of the 
structure function on a subset of our results. 

We also tested the influence of the AD power-law slope a on our 
results. As an alternative approach to keeping a fixed over the 
whole time and flux range, we allowed for a varying a. Since 
the continuum emission fr om the AD is found to get harder as 
the AD brightens (see e.g. iTrevese et al.l (I2001I) and references 
therein), we performed an alternative IBB fit allowing for a vary¬ 
ing power-law slope of the form amrit) = a ai ■ L^{t)l{L^(t)). 
We obtained best fit values a = 1.56 + 0.07 and ai =0.10 + 0.07, 
while all other parameters stayed nearly the same (Table |5]l. At 
this point, we note that unfortunately the power-law slope is only 
well-determined in those epochs where optical measurements 
were also available (epochs 7-9). Thus, our applied approach 
of fitting a varying a to our data might be insufficient to de¬ 
termine the true variability range of the AD power-law slope. 
Therefore, we alternatively used published empi rical relations 
between a and the AD luminosity. According to iTrevese et al.l 
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Fig. 10. Marginalized posterior probability distributions for the four parameters C 2 , and a of the IBB model fit, with the 1-dimensional 
projections shown along the diagonal, and the 2-dimensional projections in the other panels. Contours mark the 10%, 25%, 50%, 85% and 99% 
confidence intervals. 


(1200 Ih . who analyzed multiepoch data of a sample of quasars, 
there is a relation Aa - a + bt^(}ogFy) between the change of the 
spectral index and the logarithmic optical continuum flux change 
of the AD, with a = (-8.49 + 5.50) ■ 10~^ an d b = 2.55 + 0.75. 
Specifically for NGC 4151. lF^ti et al.l(ll984ll report the relation 
a - b ■ \ogFy with b ^ A. Making use of the inferred a values 
in epochs 7-9 and the flux differences of the z band signal with 
respect to the z flux in these three epochs, we applied the cited 
two models for deriving an alternative, more representative evo¬ 
lution of a over the whole time and flux range. Thus, instead of 
fitting the evolution of a to our data, which may problematic be¬ 
cause of the absence of optical data in most epochs, we now use 
those two models as input for our fit. We And our results to be 
qualitatively robust under the use of these different models (also 
see Sect. l53T i. even though values as high as a 3.3 are reached 
in the high-flux epochs. The influence of the particular choice of 
a for a subset of our results is shown in Table |7] 


5. Discussion 

5.1. Single-blackbody model 

As for our 2010 data of NGC 4151, we observe a significant 
rise of the emission and temperature of the innermost hot dust, 
following states of increased AD brightness. The hot dust 
temperature follows the AD flux with a time delay of roughly 
31 days. This indicates that the hot dust in NGC 4151 currently 
observed is simply heated up by increased AD irradiation and 
not destroyed owing to sublimation. In the case of significant 
dust destruction, one would expect an increase in the reverber¬ 
ation delay, which is not seen in our data. Obviously, the major 
part of the hot dust in the nucleus of NGC 4151 is not located 
at its current sublimation radius, but is cooler than sublimation 
temperature. 

There are strong indications that the hot circumnuclear 
dust around AGNs mainly consists of large graphite grains 
(» 0.2/um grain size), with sublimation temperatures > 1500K 
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dGaskell et al.l l2004t iKishimoto et al.l l2Q07l 1201 ih . Since our 
inferred dust temperatures do not reach 1500K (see Fig. [T] and 
Fig. [T2]|, it is to be expected that we see no dust sublimation in 
our data. 

Limitations of our method are given by our z-band sampling, 
and thus the interpolation of our input AD signal. Because, for 
example, the minimum observed in JHK around MJD 56700 is 
missed in the z-band observations (therefore in the interpolated 
AD signal), the resulting model hts the JHK fluxes around MJD 
56700 very poorly. Apart from this mismatch, the IBB model al¬ 
ready hts the data remarkably well within the errors. Only in the 
epochs 4-6 are the modeled ^f-band Huxes systematically lower 
than the actually measured JT-band peak huxes. 

5.2. Two-blackbody model 

In the 2BB model, epochs 4-6 are ht better, indicating that a 
second BB component might actually contribute signihcantly to 
the observed ^f-band Hux. However, we must point out that the 
goodness of the ht is not improved when using the 2BB model 
(;if^^^=1.89) instead of the IBB model (x^^^-1.87)). This is not 
surprising, since the weight of the few 2010 HK data points is 
low compared to the global data set. 

In our Paper I, we already argued that the hot, NIR dust around 
NGC 4151 might be substantially radially extended and better 
represented by more than one blackbody component. 

However, when including data at longer wavelengths, i.e. the 
WISE W1-W3 bands, we must see that our 2BB model has to 
be rejectecfl From Figs.fl4landfT5l it becomes apparent that the 
Huxes are ht better with the IBB model. 

This hnding agrees with the data of iBurtscher et al.l (l2009ll . 
The unresolved point source that they measure with A-band 
interferometry and which they attribute to the inner rim of the 
hot dust torus (at 0.04 - 0.05 pc a; 50 light-days) is apparently the 
same source that constitutes our observed NIR Huxes, because 
the observed Burtscher PS Huxes are consistent with the SED 
of our IBB model (Eig. [T4ll. In Eigs. [T4]and[T5] we went on to 
ov erplot the huxes of an extended Gaussian source measured 
by IBurtscher et akl (l2009l) . and the sum of both components (PS 
-I- Gaussian). The extended source is interpreted as the warm 
component of the clumpy torus located farther out, at 2.0 + 
0.4 pc » 240 light-days. Because this extended source is not 
resolved in the WISE photometry, it contributes to the WISE 
PSE Hux of the W3 band (which matches the Burtscher total 
Hux at 12.5 //m), while the W1 and W2 bands obviously show 
no signihcant warm dust contribution. A clear excess of the 
observed PS -H Gaussian Huxes over the modeled IBB fluxes 
can only be seen from the N band (A > 8/im) on. Here, the 
single blackbody approximation is no longer sufficient, but a 
second blackbody component contributes to, or even dominates, 
the measured fluxes. This is, however, not the 698 + 64K 
blackbody component from our 2BB Ht, but a blackbody of 
lower temperature (T = 285^5q) that is located farther out 
(IBurtscher et al.l[2009l) . A second blackbody component within 
a 100 light-day distance from the source (which we set as an 
upper limit in the prior for T 2 , see caption of Table |4]i is thereby 


* The WISE data were not included in the of the fit, because the 
aim was to test for an extended structure in the NIR regime up to the K 
band. Nevertheless, any resulting 2BB model would have to match the 
WISE fluxes as well. 


clearly ruled out. 


5.3. Variable time lag 

Besides a second blackbody component contributing to the NIR 
fluxes, there might be another explanation for the missed 2010 
A-band peak by our fit. The resulting lag of ti = 31.0+1.6 
days found for the IBB model Hts our global data very well. 
This holds for all data except for the A-band peak in epochs 
4-6. Here, it looks as if the actual lag might still be a bit 
higher. If we split the data set into two subsets (set 1: epoch 
1-9 (2010 - beginning of 2012 data), set 2: epochs 7-29 (2012 
- 2014 datafl), we And that the fltted hot dust reverberation 
delay is = 42.5 + 4.0 days if we include only epochs 1-9 
(see Table |5] for the other para meter values), consistent with 
iHonig & K^imot^ ( 1201 Ih and IKishimoto et akl (1201 lab . Eor 
the epoch 7-29 Ht, we get a best-flt delay of = 29.6 +1.7 
days (also see Table [S]). We excluded that this decrease in r 
is merely an effect of the improved sampling in the second 
period. The blackbody constant also changes between the two 
different Hts, from Cf= 4.63 + 0.33 to Cf = 3.52 + 0.22, 
and the initial temperature from Tq®'® = 1392 + 20K to 
7^013 _ ^ However, refers to the initial 

temperature in epoch 1 and thus has no real physical meaning 
for the second-period Ht, since only epochs 7-29 are included. 
Eigure [12] shows the resulting temperature evolution of this Ht. 
Interestingly, while all other parameters stay nearly the same, 
the variability factor also changes between the two fits, from 
^2010 _ Q gg ^ Q Q 4 jQ ^2013 = 0.81 + 0.03, which might indicate 
that the AD illumination is less efficiently reprocessed by the 
hot dust in the second period. 

As discussed in Sect. H] we tested the influence of differ¬ 
ent structure-function models and different evolutions of the 
AD power-law index a on our results. It seems particularly 
important to test the robustness of the apparent decrease in 
the reverberation lag under the different models. While the 
parameter values are very insensitive to the different structure 
function and a models for the global IBB model, which includes 
all epochs, the parameters inferred for the epoch 1-9 fit and 
epoch 7-29 fit do undergo certain changes when applying the 
different structure functions and a models. The influence of the 
various structure function models on our results are shown in 
Table |6| while the influence of a is shown in Table Q We do see 
a change in the derived absolute parameter values, especially 
the time lag seems to be sensitive to the applied variability 
model of a. Nevertheless, relative to each other, the inferred lags 
for the epoch 1-9 and epoch 7-29 fits and the decrease in the 
delay remain unchanged. While for a constant a (and different 
structure functions), the lag decreases from » 43 days to 
t2013 _ 3 Q (Jays, we infe r 37 days, 28 days using 

the model according to iTrevese et al.l (l200lh. and ~ 3 4 
days, ss 26 using the one presented by lEanti et al.1 (Il984l) . 
Thus, our results are qualitatively unaltered: we see a significant 
decrease in the reverberation delay from 2010 to 2013-2014. 

In Eig.|7| the temperature evolution resulting from the epoch 
1-9 and epoch 7-29 fits is plotted versus the temperature from 

^ We split whole data set into these two overlapping data sets because, 
while it was obvious that the delay in epochs 1-6 would be longer than 
the delay in epochs 10-29, a priori the delay of epochs 7-9 was not 
apparent, because fitting only 3 epochs is rather ambiguous. 
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Fig. 11. Marginalized posterior probability distributions for the four parameters C 2 , Tq,t, and a of the IBB model fit, for the first part of the 
data set (epochs 1-9) in green and the second part (epochs 7-29) in red. The 1-dimensional projections are shown along the diagonal, and the 
2-dimensional projections in the other panels. 


the single-epoch decomposition (Eq. B- Naively one would 
expect that due to the decreased AD radiation by 50% on aver¬ 
age from 2010 to end of 2012-2014 ~ 0.49), 

the blackbody temperature would have decreased by a factor 
0.5 0"^ a! 0.84 following L cx This temperature decrease 
due to the dimmer AD would be roughly balanced by a tem¬ 
perature increase due to a new dust location farther inside by 
a factor of ~ 1.18. Thus, one would expect the 

blackbody temperature for 2010 to be on the same average 
level as at the end of 2012-2014. Nevertheless, we observe a 
decreased average temperature in the second part of the data 
set ((7’)20i3^^y-^20io _ 0.86) in the single-epoch hts. This is 
conhrmed by the temperature curve of our multiepoch ht for 
epochs 7-29 and achieved through a lower variability factor 
with =» 0.84. Obviously, in the second part of the 

data set, the dust is heated less efficiently by the AD radiation 
estimated from the z band flux. The physical reasons for this 
could be an increase in dust grain size or a geometrical cause. 


In a model proposed by ICzernv & Hrvniewiczl (1201 ih . the dust 
clouds resulting from an AD wind only become exposed to 
the AD radiation in the process of moving farther outside, 
because dust that is still located farther inside tends to rise 
only a short distance above the disk, whereas the dust farther 
outside has a greater height. The innermost dust would then 
be heated less efficiently assuming an anisotropic radiation 
characteristic of the AD. The discussed parameter changes be¬ 
tween the two periods might point to a changed dust distribution. 

In Fig. [TTl we show the marginalized posterior probability 
distributions of the parameters C 2 ,Tq,t, and a for the epoch 
1-9 fit and the epoch 7-29 fit . Clearly, we see a significantly 
decreased reverberation lag from = 42.5 + 4.0 in 2010 
to = 29.6 + 1.7 in 2012-2014. As already mentioned in 
Sect, m the pdf of is bimodal. Interestingly, no bimodality 
is observed for Here, the pdf seems perfectly unimodal. 

Possibly, the bimodality is washed out by the higher photo¬ 
metric errors in the first period. Increasing the errors in the 
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Table 6. Global mean and errors for the parameters in the first period (epochs 1-9) and in the second period (epochs 7-29) for the IBB model. The 

errors are given by Ax = for each parameter. In particular, we show the influence of different structure function models discussed in Sect. l3.2l 
on our results. Here, DRWl refers to the damped random walk model with cr=0.15 and t=230, DRW2 to the damped random walk model with 
cr=0.03 and T=10. All results are shown for constant a. 


Parameter 

Unit 

ep. 1-9 
power law 

ep. 1-9 
DRWl 

ep. 1-9 
DRW2 

ep. 7-29 
power law 

ep. 7-29 
DRWl 

ep. 7-29 
DRW2 

Cl 

- 

0.78 + 0.01 

0.78 ± 0.01 

0.78 + 0.01 

0.78 + 0.01 

0.78 + 0.01 

0.78 + 0.01 

a 

- 

1.65 + 0.07 

1.65 + 0.08 

1.65 + 0.08 

1.63 + 0.04 

1.64 + 0.06 

1.64 + 0.05 

Cl 

10 ^'^ster 

4.63 + 0.33 

4.66 + 0.35 

4.67 + 0.36 

3.52 + 0.22 

3.55 + 0.34 

3.55 + 0.32 

To.i 

lOOOK 

1.392 + 0.020 

1.393 + 0.027 

1.389 + 0.026 

1.455 + 0.018 

1.449 + 0.017 

1.443 + 0.019 

Tl 

days 

42.5 + 4.0 

42.4 + 4.6 

42.2 + 4.4 

29.6+ 1.7 

30.1 + 2.1 

30.2 + 2.0 

Vl 

- 

0.96 + 0.04 

0.91 + 0.03 

0.91 + 0.04 

0.81 + 0.03 

0.77 + 0.03 

0.76 + 0.03 


Table 7. Global mean and errors for the parameters in the first period (epochs 1-9) and in the second period (epochs 7-29) for the IBB model. The 

errors are given by Ax = for each parameter. In particular, we show the influence of the evolution of a (see Sect. IS on the robustness of our 
results. All results are shown for the power-law structure function model. 


Parameter 

Unit 

ep.1-9 
a =const 

ep. 1-9 

a (Trev. 2001) 

ep. 1-9 
a (Fan. 1984) 

ep. 7-29 
a =const 

ep. 7-29 
a (Trev. 2001) 

ep. 7-29 
a (Fan. 1984) 

Cl 

- 

0.78 + 0.01 

0.75 + 0.01 

0.70 + 0.01 

0.78 + 0.01 

0.79 + 0.01 

0.77 + 0.01 

a 

- 

1.65 + 0.07 

- 

- 

1.63 + 0.04 

- 

- 

Cl 

lO-'^ster 

4.63 + 0.33 

4.83 + 0.43 

4.42 + 0.28 

3.52 + 0.22 

3.80 + 0.35 

3.70 + 0.21 

To,i 

lOOOK 

1.392 + 0.020 

1.405 + 0.036 

1.451 +0.018 

1.455 + 0.018 

1.463 + 0.027 

1.505 + 0.039 

Tl 

days 

42.5 + 4.0 

36.5 + 5.3 

33.9 + 4.2 

29.6+ 1.7 

28.3 + 2.8 

25.9 + 2.5 

Vl 

- 

0.96 + 0.04 

1.09 + 0.04 

1.15 + 0.05 

0.81+0.03 

1.03 + 0.02 

1.24 + 0.04 



Fig. 12. Hot dust temperature derived from the single-epoch decomposition, where T is a free parameter in each epoch, versus the temper¬ 
ature from the multiepoch fits, in which T evolves according to Eq.|4]and only the initial temperature To is a free parameter. For epochs 1-9, 
the temperature resulting from the epoch 1-9 fit is plotted, for epochs 10-29 the temperature from the epoch 7-29 fit. The minimum that is seen in 
r,„,gfe around MJD 56700 (and also in the corresponding JHK flux minima) is not covered within the sampling of the z-band data (see Figs.[^and 
O, therefore we observe a substantial discrepancy between Tsingle and Tm,,/,, around that date. 


second period to the level of the epoch 1-9 errors, however, 
only slightly washes out the bimodality but does not make it 
vanish completely. We furthermore excluded that the detected 
bimodality in epochs 1-9 in contrast to the missing bimodality 
in epochs 7-29 is merely due to the improved sampling in the 
second period. This supports our supposition that we see a 
changed dust distribution in the second period. 


The resolved bimodality in is caused by slightly 

different delays for the bands JHK (constituting the black- 
body emission), as confirmed by performing separate fits for 


each bancfl i.e. rf^ 26.1 + 4.8, = 28.9 + 2.6 and 

“ 31.4 + 2.4, pointing to a slight radial extent of the dust 
within a narrow region around Interestingly, what we 

observe is not just a smoothly extended structure, but at least 
two separate blackbodies at two discrete, nearby, but different 
radii. This feature is visible in Fig. [T3j especially in the H band. 
While the two-dimensional projections are hard to disentangle 
by the eye, in the right panel of Fig. one can nicely see 
an "inflection point" of the two different reverberation delays 
around t ^ 30 days. The pdf peaks at two different delays. 
Tinner ~ 26 days and Tauter ~ 31 days. In the K band. Tauter 


^ We fixed Ci,a- and v to the best fit values of the epoch 7-29 fit for 
these single-band fits, to avoid high degeneracies between the parame¬ 
ters. 
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To/IOODK r/days 

Fig. 13. Marginalized posterior probability distributions for the param¬ 
eters To and t of the epoch 7-29 IBB model fit, obtained by fitting 
the J, H, and K hands separately (also see text of this section). The J 
band fit is shown in blue, H band in green, and K band in red. The 1- 
dimensional projections of the pdfs are shown along the diagonal, and 
the 2-dimensional projections in the other panel. 


dominates the pdf, while Tinner is slightly visible as well. In the 
J band, the pdf is clearly dominated by the shorter lag, while 
the longer lag can only be marginally "resolved" due to a slight 
asymmetry in the pdf. The intersection point between the pdfs 
of the two different delays always occurs at r a; 30 days, for 
all three bands, and H and K even exhibit a minimum there. 
As expected, the fitted blackbody temperatures of the three 
bands decrease with wavelength, from - 1521 + 105K, 
= 1392 + 74K, and = 1351 + 102K. Due to missing 
information on the color in the single-band hts, the degeneracies 
between Tq and C 2 are extremely high, so the htted temperatures 
should not be taken too seriously. The detected 2BB structure 
within a narrow range around does not contradict the result 
of our previous 2BB model fit, which was mainly performed 
to test for a significantly radially extended dust distribution 
rather than a compact distribution. The two distinct blackbod- 
ies resolved within the epoch 7-29 fit still represent a fairly 
compact dust distribution, since they are very close, and are 
therefore not in conflict with our finding from the 2BB model fit. 

Our results seem to indicate a decreased time delay from 
2010 (MJD 55300-55400) to 2012-2013 (MJD 56000-56300). 
This stands in st r ong c ontrast to interferometric observations by 
iKishimoto et al.l (1201 3h . who measured an increased delay, from 
roughly 40-50 days in 2010 to 70-80 days in 2012 - a delay that 
is clearly not consistent with our data. 

It is possible that this apparent shift in the time lag only rep¬ 
resents the complicated and dynamic dust morphology with the 
dust being clearly extended and with clouds moving turbulently 
within this dust distribution. It seems possible that even slight 
changes in the dust distribution through motion of single clouds 
could cause noticeable changes in the measured torus response, 
e.g. by means of shielding. However, the shift in time lag could 


also mean that the dust radius has indeed decreased by inflow of 
dust from outside due to the low luminosity state of the AD and 
thereby decreased dust sublimation radius. Alternatively, BEL 
clouds could have been launched from farther inside due to de¬ 
creased AD luminosity. Dust condensation only occurs when the 
clouds have expanded to roughly three times their initial radius, 
which happens at a distance d =» 1000 ■ t/n, with do the initial 
distance of the clouds from the source (lElvis et al.ll2002h . Dust 
condensation would set in after roughly three to nine years from 
the initial launching of the clouds. Interestingly, our measured 
decrease in time lag seems to continue the decline from t ^ 60 
days observed end of 2008 ( Pott et al]l201C ) to t a; 50 - 45 days 
in 2009/2010 (IKishimoto et^.ll2009 . 2011a ). This total decrease 
in T from 2008-2012 seems to track a decline in the AD flux 
from 2003 un til 2006 with a delay of roughly five to six years 
(see Eig. 3 in IKishimoto et alJ (12013h . thus perfectly matching 
the timescale expected within the accretion wind scenario. 

An upper limit for the infall velocity of dust produced by stars, 
which is moving in from outside, is the free fall velocity v/ree = 


■\I2GM! Rdnsf onto the central black hole {Mbh ~ 4.6 ■ lO’’Msun, 
iBentz et ^l2006l) . We obtain an infall velocity of Vfree ~ 0.01c 
at the location of the hot dust torus. It would thus take > 1300 
days to decrease the dust radius from 43 to 30 days. This time 
span seems only marginally possible within our estimate. Al¬ 
though these qualitative arguments seem to favor the accretion 
wind scenario over the inflow model, we note that detailed mod¬ 
eling of the different scenarios is needed to reliably conclude on 
the dust formation mechanism in AGNs. 

Our measured 2010-2012 and 2012-2014 dust radii are inferred 
with one and the same method; however, the decrease over 
time of these measured radii will not fit the generally expected 
Rdust '^Lad relation, because it cannot be reasonably related 
to a decay in the momentary AD br i ghtne ss delayed by r (as also 
visible in Eig. 3 in iKishimoto et al.l (1201 3l) for other observations 
of NGC 4151 ). A variable time lag (that possibly traces the av¬ 
eraged AD signal of several years before, but not the momentary 
AD signal) might help explain the intrinsic scatter found in the 
radius-luminosity relation of dust around AGNs. 



X /yu.m 


Fig. 14. Resulting fluxes from the IBB model fit in the high tempera¬ 
ture epoch 5, plotted over wavelength. Overplotted are our derived NIR 
and MIR fluxes {zYJHK and Wl, W2, W3). The IBB model r oughly 
matches the MIR point source fluxes, given bv iBurtscher et al.l l l2009l) . 
Since the model was fit to our zYJHK data alone (also see Sect. [U, 
these are marked with green circles, while the other data are represented 
by triangles. 
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Fig. 15. Resulting fluxes from the 2BB model fit in the high temperature 
epoch 5, plotted over wavelength. Overplotted are our derived NIR and 
MIR fluxes (zYJHK and Wl, W2, W3). While matching the observed 
HK fluxes better than the IBB model, the 2BB model total fluxes by far 
exceed the WISE Wl, W2 measurements. A second blackhody compo¬ 
nent inside of 100 days can be clearly rejected. Again, the model was 
only fit to our zYJHK data (marked with green circles, while the other 
data are represented by triangles). 


6. Conclusions 

We presented updated results from our 2010-2014 NGC 4151 
photometric and spectroscopic data, which are part of our AGN 
hot dust reverberation project to monitor the evolution of the hot 
dust temperature and reverberation lag around AGNs. Our find¬ 
ings for NGC 4151 are; 

- Dust sublimation: If dust sublimation occurred in response 
to increased AD flux, it would happen at the inner edge of 
the dust distribution, thereby increasing the time delay. Al¬ 
though the AD brightness increased substantially within the 
observed time range, we see no signatures of any dust sub¬ 
limation traced by our data. In particular, after detected AD 
flux increases by more than 30% in 2010, and by roughly 
100 % from March 2012 to November 2012, we observed no 
increase in the hot dust reverberation delay, thus ruling out 
significant dust destruction for the time range 2010-2014. It 
seems that the hot dust in this galaxy is currently located be¬ 
yond its sublimation radius. 

- Dust temperature: The hot dust emission, and moreover, the 
host dust temperature in NGC 4151 closely tracks the AD 
flux variations on a short-term response timescale, which is 
roughly one month on average. We measured maximum hot 
dust temperatures lower than 1500 K throughout 2010-2014. 
The large graphite dust grains that are typically assumed for 
the hot circumnuclear dust around AGNs have a sublimation 
temperature of T^ut ^ 1500/r. It thus seems perfectly consis¬ 
tent that we saw no dust sublimation in the observed period. 

- Dust distribution radius: On a long-term timescale, we mea¬ 
sured a change of the hot dust reverberation delay of »13 
days in two years. Detailed comparison of our results with 
seemingly comparable interferometric experiments reveals 
different variations in the derived dust radii. These apparent 
observational inconsistencies could imply that the real dust 
distribution in NGC 4151 is fairly complex, so that the com¬ 
parison of hot dust radii measured with interferometry and 
reverberation is not straightforward. 

- Dust distribution: In our 2010 data, we saw slight indica¬ 
tions of a second blackhody component of » TOOK, which 
could not be confirmed without broader wavelength cover¬ 
age. Here, we presented a new analysis, now including WISE 


photometry, which rules out significant thermal radiation at 

TOOK from radii smaller than 100 days. 

To sum up, we found a decreased reverberation radius of the 
hot, circumnuclear dust in NGC 4151. Dust destruction in the 
observed epoch seems highly improbable from our data, but a 
slight change of the dust morphology seems likely. While the ob¬ 
served decrease in the reverberation delay matches the timescale 
expected within the accretion wind scenario perfectly, a radius 
decrease due to the inward motion of dust from outside seems 
only marginally possible, as estimated from an upper limit for 
the infall velocity of dust produced by stars (see Sect. l53T l. From 
our analysis, new dust formation in a cooling BLR wind ap¬ 
pears to be more likely than a radius decrease due to inflow. 
We emphasize, however, that detailed modeling is indispensable 
to reliably distinguish the different dust formation scenarios in 
AGNs. In our AGN host dust reverberation project, we moni¬ 
tor roughly 25 additional bright and variable Seyfert 1 AGNs 
with the GROND camera {grizJHK bands, ESO La Silla) and in 
the optical with the A ll-Sky Automated Survey for Supernovae 
(IShappee et alJl2014l) . It will be interesting to see whether our 
results for NGC 4151 apply to other AGNs as well. 
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Appendix A: Gelman convergence diagnostics 


is an estimator of the true target variance cr^ and is calculated by 
a weighted mean of the between-chain variance Bjn 


1 ^ 

Bjn =-- V(% - x.f ; 

m-l ^ 


J=i 


(A.3) 


i.e. the variance between the m chain means xj, (where x„ denotes 
the global mean over all chains and iterations), and the pooled 
within-chain variance W: 




1 




^ \2 


m(n — 1) 


2j 


M '=1 


(A.4) 


The true target mean /r is estimated by the global sample mean 
of X, i.e. ft - X , and the term Bl{mn) in Eq. lA. II refers to the 
sampling variability of //. As W always underestimates the true 
variance (W —> cr^ for n —> oo) for any finite n, the PSRF R will 
always overestimate the true scale reduction factor and can be 
used as convergence diagnostic: 

High values of R, i.e. values significantly above 1, indicate that 
further simulations may lead to an improved inference of the 
target distribution - either because the variance estimates in the 
numerator of Ea. IA.ll can be decreased further by running more 
iterations or because W will increase with continuing iterations, 
because the chains have not yet covered the complete target dis¬ 
tribution. If R is close to 1, it can be assumed that W has nearly 
converged to the target variance and that each of the m chains of 
n iterations is sufficiently close to the target distribution. 

The above holds if the sampler is started with a sufficiently 
overdispersed starting distribution. Ot herwise R values close t o 
1 might falsely diagnose convergence (iBrooks & Gelmanll998l) . 
Therefore, it is strongly recommended to graphically inspect the 
evolution of W and V, to make certain that they are not still 
evolving. 


We che ck for convergence using the Gelman convergenc e diag¬ 
nostics (iGelmm & Rubinll992HBrooks & Gelmanll998l) . Given 
m independent chains and 2n iterations of the sampler, we can 
calculate for the second half of the iteration^ in each step and 
for each parameter estimand x the potential scale reduction fac¬ 
tor (PSRE) 


V + Bl(mn) 

W ~ W 


(A.l) 


which serves as an (over-) estimate of the true scale reduction 
factor/? = Yl(p-. Here, 


, n — 1 B 
= -W-H - 


(A.2) 


^ The DE-MC proposal scheme presented in Eq. M might at first 
glance seem to violate one of the basic assumptions fo r moni ¬ 
toring convergenc e with the R-statistic of iGelman & RubinI (Il992h : 
IBrooks & Gelmaril jl998h . namely the assumtion t hat the m individ- 
ual chains are independent of each other. However, iTer BraaS ( l20Q6ll 
demonstrates that the conditional stationary pdf of one individual chain 
does not depend on the states of the other chains and is identical for all 
chains, so that the joint stationary pdf p(X[, ...,x,„) of the m chains is 
given simply by the product p(xi) x ... x p(x„,). Thus, the states of the 
individual chains Xi, ...,x,„ are i ndependent of each other at any itera¬ 
tion s tep after the bum-in phase dTer Braakll2006l : lMengersen & Roberll 
1200 3li . i.e., after the algorithm has become independent of the start¬ 
ing distribution. Convergence of the DE-MC algorithm can therefore 
be monitored using Gelman’s convergence criterion. 

® Only the second half is evaluated in order to lower effects of the 
influence of the starting distribution, i.e., the bum-in phase is discarded. 
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